## Load objects
stayer_coef <- rio::import("est/stayer_analyses_model_coefficients_main.csv") %>%
  dplyr::filter(model_id %in% c(209, 224)) %>%
  dplyr::select(-V1)
stayer_desc <- rio::import("est/stayer_analyses_descriptives_main.csv") %>%
  dplyr::filter(model_id %in% c(209, 224))
stayer_info <- rio::import("est/stayer_analyses_model_info_main.csv") %>%
  dplyr::filter(model_id %in% c(209, 224))


## Harmonize variable names
colnames(stayer_coef) <- colnames(stayer_coef) %>%
  gsub("X.", "", .) %>%
  gsub("\\._", "\\_", .) %>%
  gsub(".cold_rent_sqm_cwu_imp1.1", ".1.cold_rent_sqm_cwu_imp1", .) %>%
  gsub(".cold_rent_sqm_cwu_imp2.1", ".1.cold_rent_sqm_cwu_imp2", .) %>%
  gsub(".cold_rent_sqm_cwu_imp3.1", ".1.cold_rent_sqm_cwu_imp3", .) %>%
  gsub(".cold_rent_sqm_cwu_imp4.1", ".1.cold_rent_sqm_cwu_imp4", .) %>%
  gsub(".cold_rent_sqm_cwu_imp5.1", ".1.cold_rent_sqm_cwu_imp5", .) %>%
  gsub(".cold_rent_sqm_umn_imp1.1", ".1.cold_rent_sqm_umn_imp1", .) %>%
  gsub(".cold_rent_sqm_umn_imp2.1", ".1.cold_rent_sqm_umn_imp2", .) %>%
  gsub(".cold_rent_sqm_umn_imp3.1", ".1.cold_rent_sqm_umn_imp3", .) %>%
  gsub(".cold_rent_sqm_umn_imp4.1", ".1.cold_rent_sqm_umn_imp4", .) %>%
  gsub(".cold_rent_sqm_umn_imp5.1", ".1.cold_rent_sqm_umn_imp5", .) %>%
  gsub("gtyp5urban...500k_imp1.1", "gtyp5urban...500k.1_imp1", .) %>%
  gsub("gtyp5urban...500k_imp2.1", "gtyp5urban...500k.1_imp2", .) %>%
  gsub("gtyp5urban...500k_imp3.1", "gtyp5urban...500k.1_imp3", .) %>%
  gsub("gtyp5urban...500k_imp4.1", "gtyp5urban...500k.1_imp4", .) %>%
  gsub("gtyp5urban...500k_imp5.1", "gtyp5urban...500k.1_imp5", .)

## Recode x_names in stayer_coef (to match column names)
stayer_coef <- stayer_coef %>%
  mutate(x_names = x_names %>%
           as.character() %>%
           ifelse(
             startsWith(., "gtyp5urban > 500k"),
             paste0(substr(., 1, nchar("gtyp5urban > 500k")),
                    ".1",
                    substr(., nchar("gtyp5urban > 500k") + 1, nchar(.))),
             .
           )) %>%
  mutate(x_names = x_names %>%
           as.character() %>%
           ifelse(endsWith(., "gtyp5urban > 500k"),
                  paste0(., ".1"),
                  .)) %>%
  mutate(
    x_names = x_names %>%
      gsub("\\(", "", .) %>%
      gsub("\\)", "", .) %>%
      gsub("<", "\\.", .) %>%
      gsub(">", "\\.", .) %>%
      gsub("=", "\\.", .) %>%
      gsub(" ", "\\.", .) %>%
      gsub("\\:", "\\.", .) %>%
      gsub("&", "\\.", .) %>%
      gsub("-", "\\.", .) %>%
      gsub("\\`", "", .)
  )

## Split stayer_coef by model_id, drop empty columns
stayer_coef <- stayer_coef %>%
  group_by(model_id) %>%
  group_split %>%
  lapply(select_if, all_na)

## Stack by imputations
stayer_est <- stayer_coef %>%
  lapply(.,
         function(mod) {
           varying <- mod %>%
             dplyr::select(contains("_imp")) %>%
             names() %>%
             split(
               .,
               mod %>%
                 dplyr::select(contains("_imp")) %>%
                 names() %>%
                 sapply(function(name)
                   substr(name, 1, nchar(name) - 5)) %>%
                 factor(
                   .,
                   levels = mod %>%
                     dplyr::select(contains("_imp")) %>%
                     names() %>%
                     sapply(function(name)
                       substr(name, 1, nchar(name) - 5)) %>%
                     unique()
                 )
             ) %>%
             lapply(function(names)
               sapply(names, function(name)
                 which(names(mod) == name)))
           mod <- mod %>%
             as.data.frame() %>%
             reshape(
               direction = "long",
               varying = varying,
               v.names = names(varying),
               timevar = "imp",
               times = 1:5
             ) %>%
             dplyr::select(-id)
           
           return(mod)
         })

## Split by model_id and imp
stayer_est <- stayer_est %>%
  lapply(.,
         function(mod) {
           mod %>%
             group_by(imp) %>%
             group_split()
         })

## Simulate and combine
stayer_imp_sim <- stayer_est %>%
  lapply(.,
         function (mod) {
           b_sim <- lapply(mod,
                           function (imp) {
                             x_names <- imp$x_names
                             b <- as.matrix(imp[, 4])
                             vcov <- as.matrix(imp[, 5:ncol(imp)])
                             vcov <- vcov[, x_names]
                             set.seed(20230329)
                             sim <- MASS::mvrnorm(10000L, b, vcov)
                             
                             return(list(x_names = x_names,
                                         sim = sim))
                           }) 
           x_names <- b_sim[[1]]$x_names
           b_sim <- lapply(b_sim, function (obj)
             obj$sim) %>%
             do.call(rbind, .) %>%
             apply(., 2, quantile, c(.5, .025, .975)) %>%
             round(3) %>%
             format(nbig = 2, nsmall = 3) %>%
             apply(., 2, function (x) {
               paste0("\\makecell{", "$", x[1], "$ \\\\ $", "[", x[2], ",", x[3], "]$", "}")
             }) %>% as.matrix()
           b_sim <- cbind.data.frame(x_names,
                                     b_sim)
         })

dict <- c(
  "cmr_arm_cwu" = "Local market rent (EUR/sqm)",
  "log_hinc_eq_cwu" = "Equiv. household income (log)",
  "prop_personal_hinc_cwu" = "Proportion personal income",
  "hh_prop_ecact_cwu" = "Proportion econ. active household members",
  "hh_mmb_cwu" = "Number of household members",
  "lm_part_Atypical_cwu" = "\\hspace{5mm}   Atypical employment",
  "lm_part_Inactive_cwu" = "\\hspace{5mm}   Economically inactive",
  "lm_part_Unemployed_cwu" = "\\hspace{5mm}   Unemployed",
  "lm_part_In_Education_cwu" = "\\hspace{5mm}   In education",
  "lm_part_Pensioners_cwu" = "\\hspace{5mm}   Retired",
  "cold_rent_sqm_cwu" = "Household rent (EUR/sqm)",
  "cmr_arm_cwu.log_hinc_eq" = "Local market rent (demeaned) $\\times$ household income",
  "log_hinc_eq.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income",
  "Intercept" = "Intercept",
  "log_hinc_eq" = "Equiv. household income (log)",
  "age" = "Age",
  "east" = "\\hspace{5mm}   East",
  "fem" = "\\hspace{5mm}   Female",
  "edu5Upper.Secondary" = "\\hspace{5mm}   Upper-secondary",
  "edu5Post.Secondary.Non.Tertiary" = "\\hspace{5mm}   Post-secondary non-tertiary",
  "edu5Higher.Vocational" = "\\hspace{5mm}   Higher vocational",
  "edu5Higher.Education" = "\\hspace{5mm}   Tertiary",
  "gtyp3suburban" = "\\hspace{5mm}   Suburban",
  "gtyp3urban" = "\\hspace{5mm}   Urban",
  "year_fac2015" = "\\hspace{5mm}   2015",
  "year_fac2016" = "\\hspace{5mm}   2016",
  "year_fac2017" = "\\hspace{5mm}   2017",
  "year_fac2018" = "\\hspace{5mm}   2018",
  "cmr_arm_umn" = "Local market rent (EUR/sqm)",
  "prop_personal_hinc_umn" = "Proportion personal income",
  "hh_prop_ecact_umn" = "Proportion econ. active household members",
  "hh_mmb_umn" = "Number of household members",
  "lm_part_Atypical_umn" = "\\hspace{5mm}   Atypical employment",
  "lm_part_Inactive_umn" = "\\hspace{5mm}   Economically inactive",
  "lm_part_Unemployed_umn" = "\\hspace{5mm}   Unemployed",
  "lm_part_In_Education_umn" = "\\hspace{5mm}   In education",
  "lm_part_Pensioners_umn" = "\\hspace{5mm}   Retired",
  "cold_rent_sqm_umn" = "Household rent (EUR/sqm)",
  "log_hinc_eq.cmr_arm_umn" = "Local market rent (mean) $\\times$ household income",
  "log_hinc_eq.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income"
)

stayer_imp_sim <- stayer_imp_sim %>%
  lapply(
    .,
    . %>%
      dplyr::mutate(
        type = dplyr::case_when(
          grepl("_cwu", x_names) ~ "within",
          grepl("_umn", x_names) | x_names == "log_hinc_eq" ~ "between",
          TRUE ~ "mixed"
        ),
        type = factor(
          type,
          levels = c("within", "mixed", "between"),
          labels = c("within", "mixed", "between"),
          ordered = TRUE
        )
      ) %>%
      dplyr::arrange(type)
  )

## Combined table (upper part)
upper_table <- stayer_imp_sim[[1]]
upper_table <- cbind(
  upper_table,
  NA_character_
)
which_rows <- match(stayer_imp_sim[[2]][, "x_names"], upper_table[, "x_names"])
upper_table[which_rows, ncol(upper_table)] <- stayer_imp_sim[[2]][, "b_sim"]


## Separate by effect type
make_dummy_table <- function(table, dummy_pos, dummy_labels) {
  if (dummy_pos[1] == 1) {
    new_table <- rbind(c(
      dummy_labels[1],
      rep("", ncol(upper_table) - 1)
    ),
    table[dummy_pos[1]:(dummy_pos[2] - 1),])
  } else {
    new_table <- rbind(
      table[1:(dummy_pos[1] - 1),],
      c(
        dummy_labels[1],
        rep("", ncol(upper_table) - 1)
      ),
      table[dummy_pos[1]:(dummy_pos[2] - 1),]
    )
  }
  dummy_pos <- c(dummy_pos, nrow(table) + 1)
  for (n in 2:length(dummy_labels)) {
    new_table <- rbind(
      new_table,
      c(
        dummy_labels[n],
        rep("", ncol(upper_table) - 1)
      ),
      table[dummy_pos[n]:(dummy_pos[n + 1] - 1), ]
    )
  }
  return(new_table)
}

row_pos <- c(
  which(upper_table[, "type"] == "within")[1],
  which(upper_table[, "type"] == "mixed")[1],
  which(upper_table[, "type"] == "between")[1]
)
row_labels <- c(
  "\\emph{Within effects}",
  "\\emph{Mixed effects}",
  "\\emph{Between effects}"
)

upper_table <- make_dummy_table(
  upper_table,
  row_pos,
  row_labels
)

upper_table <- upper_table[, c("x_names", "b_sim", "NA_character_")]
colnames(upper_table) <- c("Variables", "Renter model", "Owner model")

## Dummy sets
dummy_pos <-
  c(
    which(upper_table[, "Variables"] == "lm_part_Atypical_cwu"),
    which(upper_table[, "Variables"] == "east"),
    which(upper_table[, "Variables"] == "fem"),
    which(upper_table[, "Variables"] == "edu5Upper.Secondary"),
    which(upper_table[, "Variables"] == "gtyp3suburban"),
    which(upper_table[, "Variables"] == "year_fac2015"),
    which(upper_table[, "Variables"] == "lm_part_Atypical_umn")
  )
dummy_labels <- c(
  "\\emph{Labor market status (ref: Full-time employment)}",
  "\\emph{East/West residence (ref: West)}",
  "\\emph{Sex (ref: Male)}",
  "\\emph{Education (ref: $<$ Upper secondary)}",
  "\\emph{Locality (ref: Rural)}",
  "\\emph{Year (ref: 2014)}",
  "\\emph{Labor market status (ref: Full-time employment)}"
)

## Map variable names
upper_table[, "Variables"] <-
  sapply(upper_table[, "Variables"], function(name)
    ifelse(name %in% names(dict),
           dict[which(names(dict) == name)],
           name))

upper_table <- make_dummy_table(
  upper_table,
  dummy_pos,
  dummy_labels
)

## Combined table (lower part, num obs)
lower_table_n <- stayer_info %>%
  dplyr::select(starts_with("N_")) %>%
  dplyr::mutate_all(as.character) %>%
  t()
lower_table_n <- cbind(
  c(
    "Observations",
    "Individuals",
    "Postcode areas",
    "Counties/cities"
  ),
  lower_table_n
)

## Combined table (lower part, variance terms)
lower_table_s <- stayer_info %>%
  dplyr::select(model_id, starts_with("sigma_")) %>%
  tidyr::pivot_longer(-model_id) %>%
  dplyr::mutate(imp = dplyr::case_when(
    grepl(".*_imp1", name) ~ 1L,
    grepl(".*_imp2", name) ~ 2L,
    grepl(".*_imp3", name) ~ 3L,
    grepl(".*_imp4", name) ~ 4L,
    grepl(".*_imp5", name) ~ 5L)
  ) %>%
  dplyr::mutate(name = substr(name, 1, nchar(name) - 5L)) %>%
  dplyr::group_by (model_id, name) %>%
  dplyr::summarize(value = mean(value)) %>%
  dplyr::mutate(
    name = dplyr::case_when(
      name == "sigma_res" ~ "$\\sigma_{\\text{Observations}}$",
      name == "sigma_id" ~ "$\\sigma_{\\text{Individuals}}$",
      name == "sigma_plz" ~ "$\\sigma_{\\text{Postcode areas}}$",
      name == "sigma_kkz" ~ "$\\sigma_{\\text{Counties/cities}}$"
    ),
    value = format(round(value, 3), nbig = 2, nsmall = 3)
  ) %>%
  tidyr::pivot_wider(names_from = model_id) %>%
  as.matrix()

## Combine
combined_table <- rbind(
  as.matrix(upper_table),
  c(
    "$N$",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table) - 2)
  ),
  lower_table_n,
  c(
    "\\emph{Standard deviations of intercepts}",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table) - 2)
  ),
  lower_table_s
)

## TeX output
print(
  xtable(
    combined_table,
    align = c("l", "l", "c", "c"),
    label = "tab:reg1",
    caption = paste0(
      "Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations."
    )
  ),
  include.rownames = FALSE,
  include.colnames = TRUE,
  sanitize.text.function = identity,
  type = "latex",
  hline.after = c(-1,
                  -1,
                  0,
                  which(combined_table[, "Variables"] == "\\emph{Mixed effects}") - 1,
                  which(combined_table[, "Variables"] == "\\emph{Between effects}") - 1,
                  which(combined_table[, "Variables"] == "$N$") - 1,
                  nrow(combined_table),
                  nrow(combined_table)),
  size = "\\small",
  comment = TRUE,
  floating = FALSE,
  tabular.environment = "longtable",
  file = "tex/reg1_2024.tex"
)


## Load objects
stayer_coef <- rio::import("est/stayer_analyses_model_coefficients_main.csv") %>%
  dplyr::filter(model_id %in% c(254, 269)) %>%
  dplyr::select(-V1)
stayer_desc <- rio::import("est/stayer_analyses_descriptives_main.csv") %>%
  dplyr::filter(model_id %in% c(254, 269))
stayer_info <- rio::import("est/stayer_analyses_model_info_main.csv") %>%
  dplyr::filter(model_id %in% c(254, 269))


## Harmonize variable names
colnames(stayer_coef) <- colnames(stayer_coef) %>%
  gsub("X.", "", .) %>%
  gsub("\\._", "\\_", .) %>%
  gsub(".cold_rent_sqm_cwu_imp1.1", ".1.cold_rent_sqm_cwu_imp1", .) %>%
  gsub(".cold_rent_sqm_cwu_imp2.1", ".1.cold_rent_sqm_cwu_imp2", .) %>%
  gsub(".cold_rent_sqm_cwu_imp3.1", ".1.cold_rent_sqm_cwu_imp3", .) %>%
  gsub(".cold_rent_sqm_cwu_imp4.1", ".1.cold_rent_sqm_cwu_imp4", .) %>%
  gsub(".cold_rent_sqm_cwu_imp5.1", ".1.cold_rent_sqm_cwu_imp5", .) %>%
  gsub(".cold_rent_sqm_umn_imp1.1", ".1.cold_rent_sqm_umn_imp1", .) %>%
  gsub(".cold_rent_sqm_umn_imp2.1", ".1.cold_rent_sqm_umn_imp2", .) %>%
  gsub(".cold_rent_sqm_umn_imp3.1", ".1.cold_rent_sqm_umn_imp3", .) %>%
  gsub(".cold_rent_sqm_umn_imp4.1", ".1.cold_rent_sqm_umn_imp4", .) %>%
  gsub(".cold_rent_sqm_umn_imp5.1", ".1.cold_rent_sqm_umn_imp5", .) %>%
  gsub("gtyp5urban...500k_imp1.1", "gtyp5urban...500k.1_imp1", .) %>%
  gsub("gtyp5urban...500k_imp2.1", "gtyp5urban...500k.1_imp2", .) %>%
  gsub("gtyp5urban...500k_imp3.1", "gtyp5urban...500k.1_imp3", .) %>%
  gsub("gtyp5urban...500k_imp4.1", "gtyp5urban...500k.1_imp4", .) %>%
  gsub("gtyp5urban...500k_imp5.1", "gtyp5urban...500k.1_imp5", .)

## Recode x_names in stayer_coef (to match column names)
stayer_coef <- stayer_coef %>%
  mutate(x_names = x_names %>%
           as.character() %>%
           ifelse(
             startsWith(., "gtyp5urban > 500k"),
             paste0(substr(., 1, nchar("gtyp5urban > 500k")),
                    ".1",
                    substr(., nchar("gtyp5urban > 500k") + 1, nchar(.))),
             .
           )) %>%
  mutate(x_names = x_names %>%
           as.character() %>%
           ifelse(endsWith(., "gtyp5urban > 500k"),
                  paste0(., ".1"),
                  .)) %>%
  mutate(
    x_names = x_names %>%
      gsub("\\(", "", .) %>%
      gsub("\\)", "", .) %>%
      gsub("<", "\\.", .) %>%
      gsub(">", "\\.", .) %>%
      gsub("=", "\\.", .) %>%
      gsub(" ", "\\.", .) %>%
      gsub("\\:", "\\.", .) %>%
      gsub("&", "\\.", .) %>%
      gsub("-", "\\.", .) %>%
      gsub("\\`", "", .)
  )

## Split stayer_coef by model_id, drop empty columns
stayer_coef <- stayer_coef %>%
  group_by(model_id) %>%
  group_split %>%
  lapply(select_if, all_na)

## Stack by imputations
stayer_est <- stayer_coef %>%
  lapply(.,
         function(mod) {
           varying <- mod %>%
             dplyr::select(contains("_imp")) %>%
             names() %>%
             split(
               .,
               mod %>%
                 dplyr::select(contains("_imp")) %>%
                 names() %>%
                 sapply(function(name)
                   substr(name, 1, nchar(name) - 5)) %>%
                 factor(
                   .,
                   levels = mod %>%
                     dplyr::select(contains("_imp")) %>%
                     names() %>%
                     sapply(function(name)
                       substr(name, 1, nchar(name) - 5)) %>%
                     unique()
                 )
             ) %>%
             lapply(function(names)
               sapply(names, function(name)
                 which(names(mod) == name)))
           mod <- mod %>%
             as.data.frame() %>%
             reshape(
               direction = "long",
               varying = varying,
               v.names = names(varying),
               timevar = "imp",
               times = 1:5
             ) %>%
             dplyr::select(-id)
           
           return(mod)
         })

## Split by model_id and imp
stayer_est <- stayer_est %>%
  lapply(.,
         function(mod) {
           mod %>%
             group_by(imp) %>%
             group_split()
         })

## Simulate and combine
stayer_imp_sim <- stayer_est %>%
  lapply(.,
         function (mod) {
           b_sim <- lapply(mod,
                           function (imp) {
                             x_names <- imp$x_names
                             b <- as.matrix(imp[, 4])
                             vcov <- as.matrix(imp[, 5:ncol(imp)])
                             vcov <- vcov[, x_names]
                             set.seed(20230329)
                             sim <- MASS::mvrnorm(10000L, b, vcov)
                             
                             return(list(x_names = x_names,
                                         sim = sim))
                           }) 
           x_names <- b_sim[[1]]$x_names
           b_sim <- lapply(b_sim, function (obj)
             obj$sim) %>%
             do.call(rbind, .) %>%
             apply(., 2, quantile, c(.5, .025, .975)) %>%
             round(3) %>%
             format(nbig = 2, nsmall = 3) %>%
             apply(., 2, function (x) {
               paste0("\\makecell{", "$", x[1], "$ \\\\ $", "[", x[2], ",", x[3], "]$", "}")
             }) %>% as.matrix()
           b_sim <- cbind.data.frame(x_names,
                                     b_sim)
         })

dict <- c(
  "cmr_arm_cwu" = "Local market rent (EUR/sqm)",
  "log_hinc_eq_cwu" = "Equiv. household income (log)",
  "prop_personal_hinc_cwu" = "Proportion personal income",
  "hh_prop_ecact_cwu" = "Proportion econ. active household members",
  "hh_mmb_cwu" = "Number of household members",
  "lm_part_Atypical_cwu" = "\\hspace{5mm}   Atypical employment",
  "lm_part_Inactive_cwu" = "\\hspace{5mm}   Economically inactive",
  "lm_part_Unemployed_cwu" = "\\hspace{5mm}   Unemployed",
  "lm_part_In_Education_cwu" = "\\hspace{5mm}   In education",
  "lm_part_Pensioners_cwu" = "\\hspace{5mm}   Retired",
  "cold_rent_sqm_cwu" = "Household rent (EUR/sqm)",
  "cmr_arm_cwu.log_hinc_eq" = "Local market rent (demeaned) $\\times$ household income",
  "log_hinc_eq.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income",
  "Intercept" = "Intercept",
  "log_hinc_eq" = "Equiv. household income (log)",
  "age" = "Age",
  "east" = "\\hspace{5mm}   East",
  "fem" = "\\hspace{5mm}   Female",
  "edu5Upper.Secondary" = "\\hspace{5mm}   Upper-secondary",
  "edu5Post.Secondary.Non.Tertiary" = "\\hspace{5mm}   Post-secondary non-tertiary",
  "edu5Higher.Vocational" = "\\hspace{5mm}   Higher vocational",
  "edu5Higher.Education" = "\\hspace{5mm}   Tertiary",
  "gtyp3suburban" = "\\hspace{5mm}   Suburban",
  "gtyp3urban" = "\\hspace{5mm}   Urban",
  "year_fac2015" = "\\hspace{5mm}   2015",
  "year_fac2016" = "\\hspace{5mm}   2016",
  "year_fac2017" = "\\hspace{5mm}   2017",
  "year_fac2018" = "\\hspace{5mm}   2018",
  "cmr_arm_umn" = "Local market rent (EUR/sqm)",
  "prop_personal_hinc_umn" = "Proportion personal income",
  "hh_prop_ecact_umn" = "Proportion econ. active household members",
  "hh_mmb_umn" = "Number of household members",
  "lm_part_Atypical_umn" = "\\hspace{5mm}   Atypical employment",
  "lm_part_Inactive_umn" = "\\hspace{5mm}   Economically inactive",
  "lm_part_Unemployed_umn" = "\\hspace{5mm}   Unemployed",
  "lm_part_In_Education_umn" = "\\hspace{5mm}   In education",
  "lm_part_Pensioners_umn" = "\\hspace{5mm}   Retired",
  "cold_rent_sqm_umn" = "Household rent (EUR/sqm)",
  "log_hinc_eq.cmr_arm_umn" = "Local market rent (mean) $\\times$ household income",
  "log_hinc_eq.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income",
  "cmr_arm_cwu.gtyp3suburban" = "Market rent (demeaned)  $\\times$ suburban locality",
  "cmr_arm_cwu.gtyp3urban" = "Market rent (demeaned) $\\times$ urban locality",
  "log_hinc_eq.gtyp3suburban" = "Household income $\\times$ suburban locality",
  "log_hinc_eq.gtyp3urban" = "Household income $\\times$ urban locality",
  "gtyp3suburban.cmr_arm_umn" = "Market rent (mean)  $\\times$ suburban locality",
  "gtyp3urban.cmr_arm_umn" = "Market rent (mean) $\\times$ urban locality",
  "gtyp3suburban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ suburban locality",
  "gtyp3urban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ urban locality",
  "gtyp3suburban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ suburban locality",
  "gtyp3urban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ urban locality",
  "cmr_arm_cwu.log_hinc_eq.gtyp3suburban" = "Market rent (demeaned) $\\times$ household income $\\times$ suburban locality",
  "cmr_arm_cwu.log_hinc_eq.gtyp3urban" = "Market rent (demeaned) $\\times$ household income $\\times$ urban locality",
  "log_hinc_eq.gtyp3suburban.cmr_arm_umn" = "Market rent (mean) $\\times$ household income $\\times$ suburban locality",
  "log_hinc_eq.gtyp3urban.cmr_arm_umn" = "Market rent (mean) $\\times$ household income $\\times$ urban locality",
  "log_hinc_eq.gtyp3suburban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ suburban locality",
  "log_hinc_eq.gtyp3urban.cold_rent_sqm_cwu" = "Household rent (demeaned) $\\times$ household income $\\times$ urban locality",
  "log_hinc_eq.gtyp3suburban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ suburban locality",
  "log_hinc_eq.gtyp3urban.cold_rent_sqm_umn" = "Household rent (mean) $\\times$ household income $\\times$ urban locality" 
)

stayer_imp_sim <- stayer_imp_sim %>%
  lapply(
    .,
    . %>%
      dplyr::mutate(
        type = dplyr::case_when(
          grepl("_cwu", x_names) ~ "within",
          grepl("_umn", x_names) | x_names == "log_hinc_eq" ~ "between",
          TRUE ~ "mixed"
        ),
        type = factor(
          type,
          levels = c("within", "mixed", "between"),
          labels = c("within", "mixed", "between"),
          ordered = TRUE
        )
      ) %>%
      dplyr::arrange(type)
  )

## Combined table (upper part)
upper_table2 <- stayer_imp_sim[[1]]
upper_table2 <- cbind(
  upper_table2,
  NA_character_
)
which_rows <- match(stayer_imp_sim[[2]][, "x_names"], upper_table2[, "x_names"])
upper_table2[which_rows, ncol(upper_table2)] <- stayer_imp_sim[[2]][, "b_sim"]


## Separate by effect type
make_dummy_table <- function(table, dummy_pos, dummy_labels) {
  if (dummy_pos[1] == 1) {
    new_table <- rbind(c(
      dummy_labels[1],
      rep("", ncol(upper_table) - 1)
    ),
    table[dummy_pos[1]:(dummy_pos[2] - 1),])
  } else {
    new_table <- rbind(
      table[1:(dummy_pos[1] - 1),],
      c(
        dummy_labels[1],
        rep("", ncol(upper_table) - 1)
      ),
      table[dummy_pos[1]:(dummy_pos[2] - 1),]
    )
  }
  dummy_pos <- c(dummy_pos, nrow(table) + 1)
  for (n in 2:length(dummy_labels)) {
    new_table <- rbind(
      new_table,
      c(
        dummy_labels[n],
        rep("", ncol(upper_table) - 1)
      ),
      table[dummy_pos[n]:(dummy_pos[n + 1] - 1), ]
    )
  }
  return(new_table)
}

row_pos <- c(
  which(upper_table2[, "type"] == "within")[1],
  which(upper_table2[, "type"] == "mixed")[1],
  which(upper_table2[, "type"] == "between")[1]
)
row_labels <- c(
  "\\emph{Within effects}",
  "\\emph{Mixed effects}",
  "\\emph{Between effects}"
)

upper_table2 <- make_dummy_table(
  upper_table2,
  row_pos,
  row_labels
)

upper_table2 <- upper_table2[, c("x_names", "b_sim", "NA_character_")]
colnames(upper_table2) <- c("Variables", "Renter model", "Owner model")

## Dummy sets
dummy_pos <-
  c(
    which(upper_table2[, "Variables"] == "lm_part_Atypical_cwu"),
    which(upper_table2[, "Variables"] == "east"),
    which(upper_table2[, "Variables"] == "fem"),
    which(upper_table2[, "Variables"] == "edu5Upper.Secondary"),
    which(upper_table2[, "Variables"] == "gtyp3suburban"),
    which(upper_table2[, "Variables"] == "year_fac2015"),
    which(upper_table2[, "Variables"] == "lm_part_Atypical_umn")
  )
dummy_labels <- c(
  "\\emph{Labor market status (ref: Full-time employment)}",
  "\\emph{East/West residence (ref: West)}",
  "\\emph{Sex (ref: Male)}",
  "\\emph{Education (ref: $<$ Upper secondary)}",
  "\\emph{Locality (ref: Rural)}",
  "\\emph{Year (ref: 2014)}",
  "\\emph{Labor market status (ref: Full-time employment)}"
)
dummy_labels <- dummy_labels[order(dummy_pos)]
dummy_pos <- dummy_pos[order(dummy_pos)]

## Map variable names
upper_table2[, "Variables"] <-
  sapply(upper_table2[, "Variables"], function(name)
    ifelse(name %in% names(dict),
           dict[which(names(dict) == name)],
           name))

upper_table2 <- make_dummy_table(
  upper_table2,
  dummy_pos,
  dummy_labels
)

## Combined table (lower part, num obs)
lower_table2_n <- stayer_info %>%
  dplyr::select(starts_with("N_")) %>%
  dplyr::mutate_all(as.character) %>%
  t()
lower_table2_n <- cbind(
  c(
    "Observations",
    "Individuals",
    "Postcode areas",
    "Counties/cities"
  ),
  lower_table2_n
)

## Combined table (lower part, variance terms)
lower_table2_s <- stayer_info %>%
  dplyr::select(model_id, starts_with("sigma_")) %>%
  tidyr::pivot_longer(-model_id) %>%
  dplyr::mutate(imp = dplyr::case_when(
    grepl(".*_imp1", name) ~ 1L,
    grepl(".*_imp2", name) ~ 2L,
    grepl(".*_imp3", name) ~ 3L,
    grepl(".*_imp4", name) ~ 4L,
    grepl(".*_imp5", name) ~ 5L)
  ) %>%
  dplyr::mutate(name = substr(name, 1, nchar(name) - 5L)) %>%
  dplyr::group_by (model_id, name) %>%
  dplyr::summarize(value = mean(value)) %>%
  dplyr::mutate(
    name = dplyr::case_when(
      name == "sigma_res" ~ "$\\sigma_{\\text{Observations}}$",
      name == "sigma_id" ~ "$\\sigma_{\\text{Individuals}}$",
      name == "sigma_plz" ~ "$\\sigma_{\\text{Postcode areas}}$",
      name == "sigma_kkz" ~ "$\\sigma_{\\text{Counties/cities}}$"
    ),
    value = format(round(value, 3), nbig = 2, nsmall = 3)
  ) %>%
  tidyr::pivot_wider(names_from = model_id) %>%
  as.matrix()

## Combine
combined_table2 <- rbind(
  as.matrix(upper_table2),
  c(
    "$N$",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table2) - 2)
  ),
  lower_table2_n,
  c(
    "\\emph{Standard deviations of intercepts}",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table2) - 2)
  ),
  lower_table2_s
)

## TeX output
print(
  xtable(
    combined_table,
    align = c("l", "l", "c", "c"),
    label = "tab:reg2",
    caption = paste0(
      "Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations."
    )
  ),
  include.rownames = FALSE,
  include.colnames = TRUE,
  sanitize.text.function = identity,
  type = "latex",
  hline.after = c(-1,
                  -1,
                  0,
                  which(combined_table[, "Variables"] == "\\emph{Mixed effects}") - 1,
                  which(combined_table[, "Variables"] == "\\emph{Between effects}") - 1,
                  which(combined_table[, "Variables"] == "$N$") - 1,
                  nrow(combined_table),
                  nrow(combined_table)),
  size = "\\small",
  comment = TRUE,
  floating = FALSE,
  tabular.environment = "longtable",
  file = "tex/reg2_2024.tex"
)

## Reduced table
rows <- c(
  "Local market rent (EUR/sqm)",
  "Local market rent (demeaned) $\\times$ household income",
  "Local market rent (demeaned) $\\times$ household income",
  "Market rent (demeaned)  $\\times$ suburban locality",
  "Market rent (demeaned) $\\times$ urban locality",
  "Market rent (demeaned) $\\times$ household income $\\times$ suburban locality",
  "Market rent (demeaned) $\\times$ household income $\\times$ urban locality"
)
upper_table3 <- upper_table2 %>%
  dplyr::slice(1:which(Variables == "\\emph{Mixed effects}")) %>%
  dplyr::filter(Variables %in% rows) %>%
  dplyr::left_join(upper_table %>%
                     dplyr::slice(1:which(Variables == "\\emph{Mixed effects}")) %>%
                     dplyr::filter(Variables %in% rows),
                   by = "Variables") %>%
  dplyr::mutate(
    Variables = ifelse(
      Variables == "Local market rent (EUR/sqm)",
      "Local market rent (demeaned)",
      Variables
    )
  ) %>%
  dplyr::select(
    Variables,
    `Renter model.y`,
    `Owner model.y`,
    `Renter model.x`,
    `Owner model.x`
  ) %>%
  dplyr::rename(
    `\\makecell{Renter model \\\\ (Overall)}` = `Renter model.y`,
    `\\makecell{ Owner model \\\\ (Overall)}` = `Owner model.y`,
    `\\makecell{Renter model \\\\ (Locality-specific)}` = `Renter model.x`,
    `\\makecell{ Owner model \\\\ (Locality-specific)}` = `Owner model.x`
  )

combined_table3 <- rbind(
  as.matrix(upper_table3),
  c(
    "$N$",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table3) - 2)
  ),
  cbind(lower_table_n, lower_table2_n[, 2:3]),
  c(
    "\\emph{Standard deviations of intercepts}",
    "\\makecell{ \\hspace{1mm} \\\\ \\hspace{1mm}}",
    rep("", ncol(upper_table3) - 2)
  ),
  cbind(lower_table_s, lower_table2_s[, 2:3])
)

## TeX output
print(
  xtable(
    combined_table3,
    align = c("l", "l", "c", "c", "c", "c"),
    label = "tab:reg3",
    caption = paste0(
      "Short regression table, selectively reporting constitutive terms of the marginal effects presented in Figs. \\ref{fig:afd_inc} and \\ref{fig:afd_region}. Coefficients and simulation-based 95\\% confidence intervals from hierarchical linear within-between models, estimated across $M=5$ imputations. For full regression output, see Tables \\ref{tab:reg1} and \\ref{tab:reg2} in the Online Appendix."
    )
  ),
  include.rownames = FALSE,
  include.colnames = TRUE,
  sanitize.text.function = identity,
  type = "latex",
  hline.after = c(-1,
                  -1,
                  0,
                  which(combined_table3[, "Variables"] == "$N$") - 1,
                  nrow(combined_table3),
                  nrow(combined_table3)),
  size = "\\tiny",
  comment = TRUE,
  floating = FALSE,
  tabular.environment = "longtable",
  file = "tex/reg3_2024.tex"
)
